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Abstract 


The  initiation  and  early  growth  of  tpots  in  channel  and 
boundary  layer  flows  it  simulated  uting  a  three  dimentional 


spectral  code 


The  timulated  spott  show  significant 


agreement  with  available  experimental  data  for  such 
quantities  at  growth  rates  and  spreading  angles. 
Disturbances  are  introduced  into  the  center  and  edge  of  the 
developing  channel  spott  to  investigate  the  relative 
sensitivity  of  spots. 
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INTRODUCTION 


Emnons  was  the  firat  to  observe  turbulent  spots  in  a 
laminar  flow  undergoing  transition  to  a  turbulent  flow. 
Since  then  a  large  number  of  investigators  have  recognized 
the  importance  of  spots  in  the  study  of  both  transition  and 
turbulence.  Naturally  occurring  spots  are  initiated  by  flow 
disturbances  like  noise.  In  the  laboratory,  spots  may  be 
artificially  initiated  with  electric  sparks  or  by  injecting 
a  jet  of  fluid.  In  a  numerical  simulation  of  spots, 
controlled  disturbances  may  be  imposed  on  a  solution  of  the 
Navi er -Stokes  equations. 

2 

Soon  after  Emmons'  discovery.  Elder  noted  that  spots 
tend  to  grow  independently  of  one  another,  even  when  they 


overlap.  Gaster' 


studied  the  linear  growth  of  small 


amplitude  disturbances  into  a  wave  packet  using  both 

laboratory  experiments  and  theoretical  analysis.  His 

theoretical  predictions  have  been  confirmed  by  laboratory 

observations  so  long  as  nonlinear  effects  are  not  important.  / 
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Wygnan8ki ,  Sokolov,  and  Friedman  conducted  an  experimental 
study  of  spots  in  a  boundary  layer.  Using  conditional 

sampling  techniques,  they  mapped  out  the  geometry  and  growth - — 

■n  For 

rates  of  a  spot  as  it  develops  in  a  boundary  layer.  Gad-el- ui 
Hak,  et  a  1 . ^  conducted  flow  visualization  experiments  on 
boundary  layer  spots  by  injecting  dye  upstream  of  the  spot 
initiation.  They  divided  the  spot  into  five  regions(see  . 
Figure  l).  Region  I  within  the  spot  overhangs  region  II,  Uv 

■  .  t.  :  „ 
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the  laminar  boundary  layer  below  the  head  of  the  spot. 
Region  III  appears  similar  to  a  turbulent  boundary  layer. 
In  regions  IV  and  V  the  flow  returns  to  a  "calm”  state.  The 
photograph  in  Figure  2  illustrates  the  characteristic 
arrowhead  shape  of  a  boundary  layer  spot  in  streamwise- 
spanwise  projection.  This  photograph  was  obtained  by 
illuminating  dye  lines  with  a  sheet  of  light  very  close  to 
the  wall. 

The  first  detailed  research  directed  toward 
investigating  the  characteristics  of  spots  in  a  channel  was 


conducted  by  Carlson,  et  al 


Using  mica  flakes  to 


visualize  the  flow  (Figure  3),  they  observed  that  a  channel 
spot  also  has  the  characteristic  arrowhead  shape.  They 
identified  (see  Figure  4)  several  features  present  in 
channel  spots.  The  spreading  half-angle  (l)  was  about  8 
degrees.  The  leading  edges  met  at  a  sharp  point  and  were 
preceded  by  oblique  waves(7).  The  center  of  the  spot  (4) 
contained  small  scale  turbulence.  Streaks(3)  trailed  from 
region  4 . 

The  purpose  of  the  present  study  is  to  use  direct 
numerical  solution  of  the  Navi e r - S t oke s  equation  to  identify 
details  of  the  internal  structure  of  spots,  as  well  as  to 
map  out  spot  dimensions  and  growth  rates.  Comparison  of  our 
results  for  growth  rates  of  the  large-scale  spot  dimensions 
with  those  seen  experimentally  verifies  that  the  essential 
growth  mechanisms  of  spots  is  captured  by  our  numerical 
expe  r iment  »  . 


One 

previous  study  of 

nume  r  i  ca  1 

spots 

should 

be 

ment ioned 

.  Leonard 

used 

discrete 

vortex 

methods 

to 

a imu  late 

suae  r  i  ca  1  ly 

the  early  growth 

of  a 

spot 

in  a 

boundary 

layer.  As  with  the 

present  computations 

,  the  spots 

computed 

by  Leonard 

are 

typical  ly 

less  mature 

than 

experimentally  observed 

spots 

• 
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2.  CXMPUTAT I ONAL  GEOMETRIES  AND  NUMERICAL  METHODS 

The  computational  domain  that  we  use  to  simulate 
channel  flow  spots  is  as  follows.  In  our  simulations  of 
channel  flow  spots,  the  flow  is  represented  by  128x64x33 
Fourier  and  Chebyshev  modes  in  the  x  ( s t reamwi s e ) ,  y 
(spanwise),  and  z  (normal)  directions,  respectively  (see 
Figure  5).  The  flow  satisfies  periodic  boundary  conditions 
in  x  and  y  and  no-slip  (rigid)  boundary  conditions  at  the 
walls  (z«±l).  The  computational  box  is  nond imens i ona 1 i zed 
by  the  channel  half-width;  in  the  runs  presented  below,  the 
physical  box  size  is  20x5x2.  With  128x64  resolution  in  x 
and  y,  the  resultant  node  spacing  (in  physical  space)  of  the 
spectral  collocation  points  is  Ax-0.16  and  Ay-0.08. 

For  our  boundary  layer  spot  calculations,  the  flow  is 
represented  using  64  Fourier  modes  in  x  and  y,  with  Ax- 2  and 
^y-l  (see  Figure  6).  In  the  z  direction,  the  33  collocation 
points  are  obtained  by  an  algebraic  mapping  of  the  interval 
[-1,1]  to  [O.oo]  with  half  the  collocation  points  located  in 
the  region  0<z<5.  The  computational  box  is 
nond  imens  i  ona  1  i  zed  by  the  boundary  layer  thickness  T]-\/vx~  /Uoo 
at  some  representative  x-location  xq .  The  periodic  boundary 
conditions  used  in  the  streamwise  direction  are  only 
approximate.  They  are  justified  because  the  increase  in 
boundary  layer  thickness  through  the  computational  domain  is 
only  67o  (see  also  Ba  1  a  sub  r  aman  i  an  ,  et  al  ).  While  inflow- 
outflow  boundary  conditions  are,  in  principle,  more 
realistic  than  periodic  boundary  conditions,  they  are  more 


r.-.v.v.v.  .".V.V, 
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wasteful  of  spatial  resolution,  which  is  the  limiting  factor 
in  the  present  calculations. 

The  Navier-Stok.es  equations  are  solved  in  rotational 

form, 

-  v  x  u  -  V(77)  +  1/Re  V2(v) 
a  t 

V-  (v)«0 

where  u>-Vxv  is  the  vorticity  and  J7-P+V  /2.  The  velocities 
are  normalized  with  respect  to  the  centerline  velocity  in 
the  channel  and  the  free  stream  velocity  in  the  boundary 
layer . 

The  spectral  method  of  Orszag  and  Patera10  is  used  in 
both  the  channel  and  boundary  layer  calculations.  For  the 
boundary  layer,  the  scheme  is  modified  by  mapping  the 
Chebyshev  collocation  points  of  the  channel  to  the  desired 
locations  in  the  boundary  layer.  A  mapping  function 

z*  -  f ( z) 

is  chosen.  'When  taking  derivatives  in  the  z-direction 
(e.g.,  in  calculating  the  vorticity)  the  Chebyshev 
differentiation  in  z*  is  followed  by  multiplication  by 
f  ’  (  z  )  : 

d d_f  d_ 

dz  ”  dz  dz* 


The  boundary  condition  at  infinity  (v  -  1)  is  implemented 
by  recalculating  the  (<  ,0)  Fourier  mode  (the  mean  flow  in  x 
and  y)  in  the  viscous  step.  Symmetry  is  not  imposed,  but 
the  spots  develop  symmetrically  when  symmetric  initial 
conditions  are  used. 

The  disturbance  is  initiated  by  applying  a  body  force 
to  a  packet  of  fluid,  producing  a  small  jet  normal  to  and 
away  from  the  wall.  The  form  of  the  disturbance  is  Gaussian 
in  x  and  y  and  continuous  in  time. 


F  -  G(t) 


l-r2/2 c2] 


where  G(t)  is  a  ramp  function.  The  size  of  the  jets  are 
indicated  in  Table  I. 


table  I 


Channe 1 


a  0.16  0.7 

Location  0.1  <  z  <  0.2  0.05  <  z  <  1.5 


Peak  normal 


Velocity 


0.09 


0.035 


We  impose  the  following  boundary  conditions  on  the  flow 
through  the  channel:  the  velocity  at  the  walls  is  zero,  and 
the  flow  is  periodic  at  the  inflow/outflow  and  cross-stream 


boundaries.  The  Reynolds  number  for  the  channel  runs  is 
6000  based  on  the  channel  half-width.  The  Reynolds  number 
for  the  boundary  layer  simulations  is  1000  based  on  the 
boundary  layer  thickness  corresponding  to  r?-1.0. 

3.  SPOTS  IN  CHANNEL  FLCKVS 

In  Figure  7,  we  plot  contours  of  the  maximum  (in  y)  of 

the  absolute  value  of  the  normal  velocity.  Max  Iv  I  for  a 

y  z 

channel  spot  at  Re-6000.  The  contour  plots  we  present  for 

channel  spots  encompass  the  entire  20x5x2  computational 

domain;  their  dimensions  are  not  to  scale.  Except  where 

noted,  the  contours  are  at  1%  intervals  of  Max  Iv  I,  where  q 

q  z 

is  the  coordinate  normal  to  the  plotting  plane.  With  this 
projection  of  the  spot  onto  a  plane  we  view  the  data  from 
the  experimentalist’s  perspective  (with  the  line  of  sight 
extending  all  the  way  through  the  channel).  At  time  t-1, 
the  initial  disturbance  has  convected  downstream  and  has 
become  slightly  distorted.  The  initial  peak  velocity  of 
0.09  has  decreased  to  0.038  due  to  viscous  diffusion.  By  a 
time  of  t-3,  the  velocity  has  increased  to  0.05  due  to 
instability  in  the  flow.  The  disturbance  is  elongated  in  x 
as  well  as  convected  downstream.  In  Figure  8  we  see  the 
disturbance  develop  most  of  the  features  characteristic  of  a 
spot.  The  front  of  the  disturbance  moves  away  from  the 
wall.  The  disturbance  grows  in  all  directions  and  the 
"arrowhead"  shape  becomes  apparent.  The  peak  normal 


velocity  increases  from  69b  at  t-12  to  99b  at  t-18. 


In  Figure  9  we  show  the  development  of  the  boundaries 
with  an  isometric  view.  Enclosed  within  the  surface  is 
fluid  whose  x  velocity  differs  from  the  Poiseuille  profile 
by  more  than  2%>. 

The  results  plotted  in  Figure  10  are  Max  Iv  I  and 

y  z 

Max^lv^l.  At  t-30,  the  spot  has  fully  extended  through  the 
channel  with  a  peak  normal  velocity  of  13%.  The  initial 
disturbance  on  the  bottom  wall  has  induced  a  new  disturbance 
at  the  top  wall.  This  second,  smaller  spot  has  a  peak 
velocity  that  occurs  at  a  distance  of  approximately  0.25 
(1/8  channel  width)  away  from  the  top  wall.  By  t-30,  the 
two  spo's  have  joined  to  produce  a  disturbance  that  fills 
the  span  of  the  channel. 

In  Figure  11  we  show  the  distortion  of  the  Poiseuille 
profile  at  the  spot  center.  The  velocity  at  the  edge  is 

essentially  unchanged  from  that  of  the  original  Poiseuille 
flow,  while  at  the  spot  center  there  is  a  velocity  defect  of 
0. 1-0.2.  At  the  bottom  wall,  the  shear  has  increased  by  a 
factor  of  3. 

In  Table  II  and  Figure  12,  we  show  how  the  spot 

geometry  changes  in  time.  Although  there  are  significant 

differences  between  conditions  generating  our  numerical  spot 

and  those  generating  the  spots  studied  experimentally,  a 

comparison  of  numerical  and  laboratory  features  is 

8 

instructive.  Carlson  et  al  generated  spots  in  a  laboratory 
channel  flow  at  Re-1000,  while  we  used  Re-6000  in  our 
calculations.  Most  of  the  experimental  data  were  taken  more 


10 


than  50-100  channel  widths  downstream  of  the  initial 
disturbance,  while  we  have  been  able  to  follow  the  spot  for 
only  10  channel  widths.  Further  development  of  the  channel 
spot  would  require  a  larger  computational  domain.  The 
growth  rate  of  the  width  and  length  of  the  numerical  spot 
becomes  constant  at  t-15  and  remains  so  until  the  spot  fills 
the  domain  at  t-32.  This  steady  growth  rate  is  slightly 
higher  than  that  observed  experimentally  in  both  the  lateral 
and  longitudinal  directions.  This  discrepancy  can  be  due  to 
the  difference  in  Reynolds  numbers  or  to  the  lack  of 
maturity  of  our  computed  spots  compared  to  those  studied 
experimentally.  We  have  not  observed  in  our  data  any 
significant  evidence  of  the  leading  To  1 lmi  en-Schl i c t i ng 
waves  that  were  observed  experimentally.  Again,  we  believe 
that  the  absence  of  these  waves  is  due  to  the  lack  of 
maturity  of  our  computed  spots. 

Table  II  Channel  Spots 

Experiment”!  I  Computational 


Velocity  of  Front  0.6 

Rear  0.34 


I  0.85 

I  0.25 


8° 


10° 


V»,a,.Va’,.vsV'i 


Spreading  Half-Angle 
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A  further  numerical  calculation  was  done  to  compare  the 
stability  characteristics  of  the  spot  at  its  edge  and 
center.  The  velocity  field  of  a  spot  at  t«20  is  used  as  the 
initial  condition  for  three  runs.  The  first  run  consists  of 
the  restarting  the  original  spot  calculation  and  allowing 
the  spot  to  continue  development  undisturbed  to  a  time  of 
24.  For  the  second  run,  a  disturbance  is  applied  at  t-20  to 
the  original  spot  at  its  center.  This  disturbance  is  of  the 
same  spatial  and  temporal  extent  as  the  original  disturbance 
that  initiated  the  spot,  but  the  magnitude  is  1/10  that  of 
the  original.  The  difference  between  the  resulting  velocity 
fields,  c(x,t)  -  I  vzl  -vz2  ••  i«  *  measure  of  the  effect  of 
the  disturbance.  By  t»24,  e(x,t)  has  exceeded  1%  in  the 
central  2/3  of  the  spot  (Figure  13).  The  third  run  is 
identical  to  the  second,  but  with  the  disturbance  applied  at 
the  spot  edge,  rather  than  at  the  center.  At  t-24  the 
disturbance  had  propagated  through  most  of  the  spot  (see 
Figure  14),  and  had  a  peak  magnitude  of  about  4%.  as  opposed 
to  the  1.5%  peak  from  the  second  run. 

From  these  results,  we  conclude  that  channel  spots  are 
less  stable  at  their  edges  than  at  their  centers.  This 
observation  suggests  that  spots  grow  by  destabilization  of 
neighboring  fluid,  rather  than  simply  engulfing  laminar 
fluid . 


4.  SPOTS  IN  BOUNDARY  LAYER  FLCKVS 


In  Figure*  15  through  17  we  *how  the  development  of  * 

boundary  layer  epot  at  Re-1000  up  to  t-90.  The  contour  plots 

we  present  for  boundary  layer  spots  encompass  the  entire 

128x64  computational  domain  in  x  and  y  and  are  truncated  at 

z-22.  Again,  except  where  noted,  the  contours  are  at  1% 

intervals  of  Max  Iv  I,  where  q  is  the  direction  normal  to 

q  z 

the  plotting  plane.  Figure  IS  shows  the  streamwise  and 
spanwise  development  of  the  spot  from  the  initial 
disturbance.  At  t-90,  the  spot  has  begun  to  develop  the 
characteristic  arrowhead  shape,  which  is  more  apparent  in 
the  second  (29b)  velocity  contour.  Figure  16  shows  the 
development  of  the  triangular  shape  and  the  overhang  in  the 
spanwise  direction.  Figure  17  shows  the  overhang  develop  in 
the  leading  edge . 

Table  III  Boundary  Layer  Spots 


Expe  r iment  a  1 

1  Computational 

Velocity 

of  Front 

0.9 

1  0.85 

Rear 

0.5 

1  0.3 

1 

Spread ing 

Half -Ang 1 e 

10° 

1 

1  12° 

The  growth  and  development  of  the  spot  in  a  boundary 
layer  is  compared  with  the  experimental  findings  of 
Wygnanski,  et  al^  in  Table  III.  The  growth  rate  of  the  spot 


in  the  itretmviie  and  apanwiae  directiona  is  in  relatively 
cloae  agreement  with  the  experimental  data.  This  auggesta 
that  the  growth  mechaniama  in  a  boundary  layer  apot  have 
been  accurately  captured  in  thia  aimulation. 


Figure 

18  shows 

cross 

sections  of  the 

spot 

at  t-90 . 

Here  we  plot 

contours 

of 

the 

local  values 

of 

vz  at  y' 

y  *  -0.5, 

■’center 

2.5,  and 

4.5, 

in 

Figures  18a, 

18b, 

and  18c, 

reapect ively .  Intervala  are  at  1%  and  daahed  contour  lines 
represent  negative  z  velocities.  The  velocities  are  highest 
in  the  plane  closest  to  the  center  of  the  apot  (see  Figure 
18a).  Away  from  the  spot  centerline  the  velocities  and  the 
spot  height  decrease.  The  front  of  the  spot  haa  an  overhang 
of  a  distance  of  10-20  in  as  has  been  observed 
experimentally.  The  flow  is  dominated  by  eddies  with  length 
scales  of  approximately  10  in  x  and  5  in  y.  These  length 
scales  differ  from  those  of  unstable  modes  of  the  Orr- 
Sonxnerfeld  equations,  which  predicts  linear  instability  for 
much  longer  wavelengths,  30<A^<85. 

In  order  to  explore  the  later  time  evolution  of 
boundary  layer  spots,  it  will  be  necessary  to  use  higher 
resolution  simulations,  which  we  hope  to  perform  in  the 
future  . 

CONCLUSIONS 


It  has  been  shown  that  spots  can  be  generated  by 
numerical  solution  of  the  Navier-Stok.es  equations.  The  fact 


that  our  results  for  the  growth  rates  of  the  large-scale 
spot  dimensions  are  relatively  close  to  those  seen 
experimentally  suggests  that  the  essential  growth  mechanisms 
of  spots  have  been  captured  by  our  numerical  experiments. 
These  simulated  spots  are  less  mature  than  typical 
experimental  spots,  but  their  behavior  appears  to 
approximate  that  in  a  fully  developed  spot. 

The  spots  generated  were  not  dominated  by  two 
dimensional  Tol lmi en-Schl i c t ing  waves.  This  suggests  that 
the  growth  in  spots  is  not  linear  growth  of  two  dimensional 
Tollmien-Schlicting  waves.  Moreover,  the  perturbation 
velocities  seen  were  about  0.1;  perturbations  this  large 
would  make  the  results  of  linear  theory  inapplicable  and 
suggest  domination  of  nonlinear  effects.  This  does  not  rule 
out  the  importance  of  Tollmien-Schlicting  waves  in  the 
amplification  of  small  disturbances  which  may  develop  into 
spots  or  as  a  driving  mechanism  for  some  secondary 
instability  in  spots. 
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FIG.  1  Schematic  of  an  experimental  boundary  layer  spot  cut 
through  the  center  ( f torn  Gad* e  1  -Halt  et  al.^). 

FIG.  2  Vi aua 1 i za t i on  of  an  experimental  boundary  layer  spot 
using  fluorescent  dye  and  a  sheet  of  laser  light  at  the 
wall;  Re^-SxlO^  (  f  tom  Gad  -  e  1  -Halt  et  al.^). 

FIG.  3  Visualization  of  an  experimental  channel  spot  using 

g 

mica  platelets  (from  Carlson  et  al.  ). 

FIG.  4  Channel  spot  schematic:  (l)  spreading  half  angle;  (2) 
trailing  streaks;  (3)  region  of  small-scale  turbulence  (4) 

g 

oblique  To  1 lmi en-Schl i c t ing  waves  (from  Carlson  et  al.  ). 
FIG.  5  Channel  geometry  and  nomenclature.  Channel  is  20x5x2 
in  the  x,y,  and  z  directions,  with  128x64  Fourier  modes  in  x 
and  y  and  33  Chebyshev  modes  in  z. 

FIG.  6  Boundary  layer  geometry  and  nomenclature.  Boundary 

layer  computational  domain  is  128x64  in  the  x  and  y 

directions,  with  64x64  Fourier  modes  in  x  and  y  and  33 

Chebyshev  modes  mapped  in  the  normal(z)  direction. 

FIG.  7  Early-time  evolution  of  channel  spot.  Max  Iv  I 

y  z 

contours  are  plotted  at  1%  intervals. 

FIG.  8  Channel  spot  at  intermediate  times.  Max^lv^l 
contours  in  a)  and  b);  Max  Iv  I  contours  in  c)  and  d). 

y  * 

FIG.  9  Surfaces  of  2?o  x-velocity  perturbations  in  developing 
channe 1  spot. 

FIG.  10  Channel  spot  at  t-30.  Max^lv^l  contours  in  a); 
Max  Iv  I  contours  in  b). 

y  * 

FIG.  11  Mean  velocity  profiles  at  center  (solid)  and  edge 
(broken)  of  spot. 


FIG.  12  Location  in  z  of  front,  center,  and  rear  of  channel 


spot  vs.  time,  where  spot  is  defined  as  region  where 
lv  l>2%.  For  t  larger  than  30,  the  spot  length  reaches  the 
periodicity  length  of  the  computational  domain,  so  the  spot 
ceases  to  grow  in  the  streamwise  direction. 

FIG.  13  Perturbation  velocity,  c(x,t),  contours  at  t-22  and 
t-24  for  channel  spot  perturbed  at  its  center  at  t-20. 

FIG.  14  Perturbation  velocity,  e(x,t),  contours  at  t-22  and 
t-24  for  channel  spot  perturbed  at  its  edge  at  t-20. 

FIG.  15  Development  of  boundary  layer  spot.  Max^lv^l 
contours  are  plotted  at  1%  intervals.  a)t-30;  b)t-60;  c)t-90 


FIG.  16  Development  of  boundary  layer  spot.  MaXylv^l 


contours  are  plotted  at  1%  intervals.  a)t-30;  b)t-60;  c)t-90 


Early-Time  Spot  Evolution  at  R 


Figure  7 
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Figure  10 
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MEAN  VELOCITY  PROFILES 
AT  CENTER  AND  EDGE  OF  SPOT 
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Numerical  Solution  of  Incompressible  Flows 
by  a  Marching  Multigrid  Nonlinear  Method 


Moshe  Rosenfeld*  and  Moshe  Israeli! 
Techruon- Israel  Institute  of  Technology,  Haifa,  Israel 


\  downstream  marching  iterative  scheme  for  the  solution  of  the  steady,  incompressible,  and  two-dimensional 
parabolized  or  thin  layer  Navier-Stokes  equations  is  described  for  a  general  curvilinear  orthogonal  coordinate 
system.  Modifications  of  the  primitive  equation  global  relaxation  sweep  procedure  result  in  an  efficient  marching 
scheme.  This  scheme  takes  full  account  of  the  reduced  order  of  the  approximate  equations  as  it  behaves  like  the 
SLOR  method  for  a  single  elliptic  equation.  Tbe  proposed  algorithm  is  essentially  Reynolds  number- independent 
and  therefore  can  be  applied  to  the  solution  of  the  incompressible  Euler  equations.  A  judicious  choice  of  a 
staggered  mesh  enables  second-order  accuracy  even  in  the  marching  direction.  The  improved  smoothing  properties 
permit  the  introduction  of  multigrid  acceleration  The  convergence  rates  are  simitar  to  those  obtained  by  the 
muliigrid  solution  of  a  single  elliptic  equation;  the  storage  is  also  comparable  as  only  the  pressure  has  to  be  stored 
on  all  levels.  Numerical  results  are  presented  for  several  boundary -layer- type  flow  problems  including  the  flow 
over  a  spheroid  at  zero  incidence. 


I.  Introduction 

CONSIDERABLE  evidence  has  accumulated  recently 
about  the  applicability  of  the  parabolized  Navier-Stokes 
(PNS)  equations  for  high  Reynolds  number  flows  with  a 
principal  flow  direction  (see  Rubin1)  The  PNS  equations  are 
obtained  bv  neglecting  the  streamwise  viscous  terms  in  the 
Navier-Stokes  equauons.  When  the  viscous  terms  in  the  cir¬ 
cumferential  direction  arc  also  neglected,  one  gets  the  thin 
layer  (TL)  approximation 

The  steady  PNS  or  TL  equations  still  have  an  elliptic  nature 
(but  of  reduced  order—see  Sec.  II  )  and  therefore  the  initial 
value  problem  in  the  downstream  marching  direction  is  not 
well  posed  -  A  well-posed  initial-boundary  value  problem  can 
be  formulated  bv  specifying,  for  example,  upstream  and  side 
conditions  for  the  velocities  and  one  downstream  condition 
for  the  pressure  Therefore  the  PNS  equations  must  be  solved 
gk'balls  and  cannot  be  solved  by  a  single-sweep  marching 
The  reduced  order  of  the  PNS  equations  can  be  exploited 
bs  constructing  an  iterative  marching  method  for  updating  the 
pressure  field  only  1  Such  a  multiple-sweep  iteration  method 
has  the  advantage  that  the  velocity  field  is  generated  during 
the  marching  process  and  only  the  pressure  field  has  to  be 
stored  from  sweep  to  sweep  — a  considerable  saving  in  storage 
results  However,  simple-minded  marching  does  not  result  in 
good  convergence  properties  and  sometimes  diverges.  For  the 
two-dimensional  case,  Israeli  and  Lin’ 4  devised  a  stable 
marching  scheme  that  behaves  like  the  successive  line  over 
relaxation  (SLOR)  method  for  the  solution  of  a  single  elliptic 
equation. 

Rubin  and  Reddy'  analyzed  certain  aspects  of  the  solution 
of  the  PNS  equations  and  used  their  procedure  to  solve 
wverui  flow  problems  in  Cartesian  and  axisymmctric  body 
titled  conformal  coordinates  They  used  in  most  cases  a  first- 
order  scheme  and  also  applied  a  simplified  one-dimensional 
mullignd  algorithm  Several  more  recent  works  are  cited  in 
Refs  fi  k 
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The  present  study  is  a  continuation  of  the  work  presented 
in  Ref.  9,  where  the  scheme  of  Israeli  and  Lin14  was  modified 
into  a  second-order,  staggered,  marching  multigrid  form.  The 
principal  aim  here  is  to  test  the  convergence  properties  of  the 
methods  in  the  case  of  viscous  nonlinear  problems  and  to 
apply  the  algorithm  to  several  flow  problems.  The  good 
smoothing  properties  of  the  Israeli  and  Lin  scheme  are  used  in 
a  multigrid  framework  in  order  to  accelerate  the  convergence 
of  the  solution  of  the  PNS  equations.  The  steady,  incom¬ 
pressible,  viscous,  and  two-dimensional  equations  are  consid¬ 
ered  in  a  general  curvilinear  coordinate  system.  The  marching 
scheme  is  implemented  using  a  stable  algorithm  which  is 
second-order  accurate  also  in  the  marching  direction.  The 
same  method  can  be  used  without  modification  for  the  incom¬ 
pressible  Euler  equations,  as  the  effect  of  the  Reynolds  num¬ 
ber  on  the  convergence  rate  is  insignificant.  In  two  dimen¬ 
sions,  the  PNS  and  the  TL  equations  are  identical,  and 
therefore  the  same  analysts  applies  to  both. 

II.  Formulation 

The  nondimensional  steady,  incompressible,  and  two- 
dimensional  PNS  (or  TL)  equations  in  a  general  curvilinear 
orthogonal  coordinate  system  (£,  tj)  are  as  follows: 

Continuity 
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4  is  approximately  aligned  with  the  mainstream  direction.  U 
and  l'  are  the  contravariant  velocity  components  in  the  4  and 
i]  directions  respectively.  P  is  the  pressure  and  Re  is  the 
Reynolds  number.  h(  and  h  are  the  Lammc  coefficients  in 
the  4  and  ij  directions.  J  is  the  Jacobian  of  the  coordinate 
transformation. 

The  two-dimensional  Navier-Stokes  equations  are  elliptic  of 
order  four  11  It  can  be  shown  that  the  PNS  or  TL  equations 
are  elliptic  of  order  two  only— like  a  single  Poisson  equation. 
This  eliipticity  is  caused  by  the  pressure  gradient  terms  via  the 
continuity  equation.  A  well-posed  problem  can  be  formulated 
by  defining  the  boundary  conditions  as  described  in  Fig.  1  for 
a  rectangular  control  volume.  The  following  conditions  may 
be  specified: 

Upstream  boundary  (AB) 


Umj,Vmr  together  with  the  pressure  Pm . ,  ,  (for  j  -  1, 
2 Thus,  the  velocities  are  indeed  solved  using  values 
from  the  upstream  while  the  pressure  uses  information  from 
the  downstream.  This  approach  was  subsequently  used  in 
Refs.  3  and  8. 

In  addition,  one  may  stagger  the  velocity  V  with  respect  to 
the  other  variables  as  shown  in  Fig.  2,  where  the  centering 
points  of  the  difference  equations  are  also  plotted.  The  dif¬ 
ferential  equations  are  approximated  by  central  second-order 
approximations.  Whenever  needed,  averaging  was  used,  as  is 
usually  done  for  staggered  grids. 

The  nonlinear  algebraic  equations  are  linearized  bv  either  a 
full  Newton- Raphson  (NR)  method  or  by  performing  only 
one  NR  iteration.  The  results  are  very  similar,  so  usually  a 
single  NR  iteration  is  used. 


v-vm  (2a) 

Solid  wall  (AD) 

(7-0  f'-O  (2b) 

Outer  boundary  (BQ 

«4m.  p  - p**  <2c) 

Downstream  boundary  (CD) 


Other  boundary  conditions  can  be  used,  but  the  same  number 
of  conditions  on  each  boundary  must  be  kept.  Subscripts  “in” 
and  “out”  refer  to  the  inner  and  outer  boundaries,  respec¬ 
tively. 

III.  Discretization 

Numerical  solutions  of  Eqs.  (1)  are  obtained  by  spreading  a 
mesh  over  the  computational  domain.  Let  us  assume  that  the 
grid  points  are  distributed  evenly  along  the  4  and  ij  coordi¬ 
nates  with  the  spacing  A  and  Atj  respectively.  When  discretiz¬ 
ing  these  equations  it  should  be  remembered  that  their  nature 
should  be  reflected  in  the  finite-difference  approximation.1-5  In 
order  to  be  consistent  with  the  boundary- layer  (parabolic) 
nature  of  the  flow,  the  axial  gradients  of  the  velocities  should 
be  computed  using  only  upstream  values,  while  the  elliptic 
nature  is  preserved  by  forward  differencing  the  axial  pressure 
gradient15  *  Consequently,  it  was  assumed  that  a  stable 
marching  scheme  must  be  of  the  first  order  in  the  marching 


IV.  The  Multigrid  Algorithm 

The  multigrid  technique  is  a  numerical  procedure  for  sub¬ 
stantially  improving  the  convergence  rate  of  iterative  methods. 
In  order  to  facilitate  comparison  with  theory,  the  accommo¬ 
dative  C-cycle  MG  algorithm  was  chosen.10  Each  MG  process 
consists  of  three  basic  parts:  relaxation,  restriction,  and  inter¬ 
polation. 

Some  of  the  elements  of  the  present  approach  were  used 
independently  in  Ref.  5.  Detailed  comparisons  cannot  be 
made  because  convergence  rates  were  not  presented  there.  In 
the  present  study  the  MG  refinement  is  applied  in  both  the  4 
and  i)  directions,  whereas  in  Ref.  5  the  computational  mesh 
was  refined  only  in  the  streamwise  direction  (one-dimensional 
MG  procedure). 

The  Relaxation  Scheme 

The  overall  convergence  rate  of  any  MG  process  is  strongly 
influenced  by  the  smoothing  properties  of  the  relaxation 
scheme.  It  can  be  shown  analytically  and  experimentally  that 
the  usual  multiple-sweep  marching  does  not  have  good  con- 


ony  two  of  (u,v,p) 


ony  two  of  (u,v,  p) 

Fig.  1  Example  of  permissible  boundary  conditions. 


direction  But  a  second-order  accuracy  can  be  achieved  by  a 
judicious  choice  of  the  placement  of  the  variables  to  be  solved 
at  each  station.  The  choice  can  be  explained  most  easily  by 
considering  a  Cartesian  coordinate  system  and  taking  V-0 
and  1  /Re  -  0  in  Eq.  (lc),  yielding 


vlm  ~pl 


(3) 


A  first-order  difference  scheme  then  becomes  (see  Fig.  2) 


u- 


U, 


m  -  1 ,  / 


'  Pm.j  Pm*\,j 


(4) 


with  the  unknowns  and  Pm,,.  An  alternative  scheme, 
first  suggested  by  Israeli'’11  is  written 


(5) 


with  the  unknowns  U„  and  P„  1(.  The  scheme  is  centered 


about  m  v  1/2  and  is  of  second  order. 

The  same  applies  to  the  full  viscous  PNS  or  TL  equations. 
At  each  marching  step  one  solves  for  all  the  velocities 


Centering  Of  Eqs.  : 

O  -  Continuity  ft  X  Momentum 
□  -  Y  Momentum 


Fig.  2  Th«  staggered  grid. 
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vergence  and  smoothing  properties  because  short  wave  errors 
are  not  efficiently  smoothed  Israeli  and  Lin1,4  showed  that 
certain  modifications  in  the  streamwise  momentum  equation 
which  vanish  upon  convergence  give  rise  to  an  iterative  scheme 
that  is  equivalent  in  the  linear  case  to  the  SLOR  method  for  a 
single  Poisson  equation  In  the  general  nonlinear  case,  the 
modified  iterative  process  is  essentially  equivalent  to  the  re¬ 
laxation  of  a  single  nonlinear  Poisson-like  equation  for  the 
pressure  The  velocities  can  be  viewed  as  auxiliary  variables 
needed  during  the  marching  since  they  have  no  “memory"  by 
themselves. 

Furthermore,  the  good  smoothing  properties  of  the  line 
relaxation  scheme  of  a  single  Poisson  equation  were  automati¬ 
cally  gained  The  problems  associated  with  the  loss  of  elliptic  - 
ity  of  the  difference  approximation  for  the  Navier-Stokes 
equations  at  high  Reynolds  number10  are  thus  avoided,  and 
no  upstream  weighting  or  artificial  viscosity  is  required.  There 
results  a  considerable  saving  in  storage  as  well  as  a  simpler 
relaxation  scheme  where  the  convergence  rate  is  essentially 
independent  of  the  Reynolds  number.  The  same  marching 
algorithm  can  thus  be  used  for  the  Euler  equation  with  the 
same  favorable  convergence  rate. 

The  extension  of  the  marching  scheme  of  Ref.  4  to  a  general 
curvilinear  orthogonal  coordinate  system  yields  a  modified 
mainstream  momentum  difference  equation  (the  other  dif¬ 
ference  equations  remaining  unchanged): 

-GmP:_l-R„  +  Sm  (6a) 

K  ,  -*/>„•_, -0  (6b) 

Rm  includes  the  velocity  terms  of  the  finite-difference  ap¬ 
proximation.  m  is  the  marching  step  number,  k  is  the  global 
iteration  index,  and  u  is  the  over-relaxation  parameter.  G„ 
and  S„,  are  given  in  Table  1.  When  the  MG  procedure  is 
applied,  the  over-relaxation  parameter  is  u  -  1. 

In  each  marching  step,  a  block  tridiagonal  system  is  solved 
for  the  vectors  Vm,  Vm,  and  each  component  of  the 

vectors  corresponds  to  a  point  along  the  ij  coordinate.  Thus, 
the  continuity  and  the  T|-momentum  difference  equations  are 
solved  exactly  in  each  step.  The  streamwise  (()  momentum 
equation  is  not  solved  exactly  since  the  pressure  Pm,  taken 
from  the  last  global  iteration,  appears  in  it. 

A  linear  von  Neumann  stability  analysis  of  the  marching 
iterative  scheme  for  the  primitive  coupled  system  of  difference 
equations  of  Eq.  (1)  in  Cartesian  coordinates  ( Z,  T )  was 
performed  by  Rubin  and  Reddy.5  They  numerically  de¬ 
termined  the  eigenvalues  of  the  matrix  of  the  relaxation  pro¬ 
cess  and  found  the  estimate  XmM  -  1 -- C2(AZ/Ym)2  for 
( A  7.  /  1  Y„,  is  the  outer  boundary  in  the  normal  T 

direction,  AZ  is  the  interval  in  the  marching  direction  Z.  The 
rate  of  convergence  is  determined  by  the  ratio  A Z/Ym.  They 
claim  that  this  conclusion  “...differs  from  that  found  for  the 


convergence  analysis  of  line  relaxation  procedures  for  Poisson 
solvers  Although  the  source  term  S,t  leads  to  the  conven¬ 
tional  relaxation  form  of  the  Laplace  operator  for  the  pres¬ 
sure,  the  coupling  with  the  velocities  alters  the  structure  of  the 
inversion  matrix  and  associated  spectral  radius." 

It  turns  out,  however,  that  the  spectral  radius  can  be 
determined  analytically  for  the  coupled  PNS  system  and  is 
independent,  in  the  linear  case,  of  the  coupling  by  the  veloci¬ 
ties.  A  stability  analysis  of  the  discretized  version  of  Eq.  (1) 
using  the  modified  form  Eq.  (6)  reveals  that  the  amplitude  of 
the  errors  in  the  velocities  can  be  related  to  the  error  in  the 
pressure  field.  On  the  other  hand,  it  was  shown5'4  that  in  the 
linearized  case  the  global  iteration  scheme  is  reducible  to  the 
SLOR  method  for  a  single  Poisson  equation.  It  follows  that  all 
the  known  results  from  the  theory  of  the  SLOR  scheme  should 
be  applicable  to  the  present  version  of  the  global  iterative 
scheme  for  the  PNS  equations.  In  particular,  we  find  that  the 
maximum  eigenvalue  X,^,,  which  determines  the  convergence 
rate,  is  given  in  the  present  case  by12 
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for  the  boundary  conditions  of  Eqs.  (2)  and  with  the  ap¬ 
propriate  scaling.  Ym  and  AT  are  the  outer  boundary  and  the 
interval  in  the  T  direction  respectively.  For  u  -  1  the  maximal 
eigenvalue  is  given  by 

^-l-H^z/rj2  (7c) 
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Numerical  experiments  verified  the  validity  of  Eq.  (7c).  For 
large  enough  A Z/ ,  ihe  convergence  is  very  rapid  and  no 
coarsening  of  the  grid  is  necessary  in  the  MG  procedure.  For 
smaller  values  of  AZ/y„,  the  MG  procedure  is  invoked  in 
order  to  accelerate  convergence.  For  very  small  values  of 
A  Z  y„,  the  MG  scheme  does  not  seem  to  be  effective,  and 
this  indicates  that  the  underlying  relaxation  scheme  breaks 
down  Since  the  linear  analysis  assures  convergence  for  all 
values  of  A Z/Y„,,  it  is  possible  that  nonlinear  effects,  which 
were  neglected  in  the  analysis,  adversely  affect  the  conver¬ 
gence  rate  and  the  smoothing  properties  of  the  relaxation. 

Restriction  and  Stonge  Requirements 

Lei  the  finite-difference  approximation  of  Eqs.  (1)  on  the 
finest  grid  V/  be  represented  as  in  Ref.  10: 

*"<*)-/■"(*)  (g) 

where  v=-i(i,tj),  W  *  -  U*,  V*.  P*]T  is  the  exact  solution 
of  the  difference  equations  and  j  is  the  number  of  the 
differential  equation  (y  -  1,2, 3). 

The  problem  is  transferred  in  the  full  approximation  stor¬ 
age  (FAS)  mode  from  the  current  level  k  to  a  coarser  level 
A  —  I  by  correcting  the  right-hand  side  of  Eq.  (8)  (see  Fig.  3): 


r  1/r‘w‘(x)] +//‘t->{F/(x)-L/‘w‘(x)} 

(9) 

w  it)  is  an  approximation  to  Wi ( x).  /‘j1  and  /*'*  are 
proper  restriction  operators  for  the  j  th  equation  and  tor  the 
dependent  variables,  respectively. 

The  term  within  the  braces  in  Eq.  (9)  is  the  residual  of  the 
j  th  equation.  For  the  present  inarching  scheme  there  is  no 
residual  in  the  continuity  and  in  the  i)-momentum  equations 
since  they  are  solved  exactly  in  each  step.  The  residual  of  the 
{-momentum  equation  results  only  from  the  streamwise  pres¬ 
sure  gradient  term,  and  its  computation  needs  only  one  sub¬ 
traction  //  1  was  chosen  to  be  a  linear  interpolation,  which 
yields  L\  1  ( /* ' 1  w‘ ( it )]  -  0  for  the  continuity  equation.  /,*  j 1 
is  computed  by  averaging  m  both  the  (  and  t)  directions. 
/‘  4 1  is  a  simple  injection.  In  summary,  Eq  (9)  takes  the 
following  form: 
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Fig.  S  Accuracy  test  of  the  pressure  and  the  normalized  skin  friction 
coefficient  at  Z  -  4  (semi-infinite  fiat  plate). 

Two  consequences  should  be  emphasized:  1)  only  two  correc¬ 
tions  have  to  be  computed  and  stored 

in  the  coarse  grids,  and  2)  oil  the  dependent  variables  must 
be  transferred  in  order  to  compute  the  corrections 
{  j  -  2,3).  Since  only  the  pressure  is  stored, 

these  corrections  must  be  computed  during  the  marching 
process. 

Special  attention  must  be  paid  to  the  restriction  of  the 
staggered  flow  variables  at  the  boundaries  in  order  to  specify 
zero  corrections  as  boundary  conditions  for  the  coarse  grids. 
In  the  present  study,  all  the  boundary  conditions  were  given  at 
the  physical  boundaries.  The  values  at  the  actual  staggered 
placement  of  each  variable  at  the  boundaries  were  computed 
by  averaging  with  inner  points. 

It  follows  that  in  addition  to  the  pressure  on  all  grids,  one 
has  to  save  one  correction  term  for  each  momentum  equation 
on  the  coarser  grids.  Assuming  N  computational  points  on  the 
finest  grid,  a  simple-minded  estimate  gives  14A/3  storage 
locations  for  the  two-dimensional  NS  multigrid  solution  and 
2  N  for  the  PNS  marching  MG  solution. 

interpolation 

Since  the  present  marching  scheme  generates  the  velocity 
field  from  the  pressure  field,  only  the  correction  to  the  pres¬ 
sure  must  be  interpolated  back  to  the  fine  grid: 

Fk  +  '-P'  +  l  +  /*♦«(/>* -/**,/>*“)  (11) 

/*"'  is  the  interpolation  operator,  fn  the  present  case  it  is  a 
linear  operator. 

The  MG  scheme  described  above  is  general.  But  so  far  only 
cases  in  Cartesian  coordinate  systems  with  equally  spaced  grid 
points  were  tested  with  the  MG  procedure. 

V.  Results 

Line*  Casa 

A  linearized  version  of  the  PNS  equations,  expressed  in  a 
Cartesian  coordinate  system,  has  been  tested  in  Ref.  9.  Some 
of  the  results  are  presented  here  for  reference  and  for  com¬ 
parison  with  the  nonlinear  solutions. 

Figure  4  compares  the  convergence  history  of  different 
relaxation  schemes  with  and  without  MG  acceleration  for  a 
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Fig.  6  Convergence  history  on  the  finest  grid  (semi-infinite  flat  plate). 


finest  grid  consisting  of  17  x  17  mesh  points.  The  horizontal 
axis  gives  the  number  of  work  units  (WU).10  A  work  unit  is 
equivalent  to  one  global  iteration  on  the  finest  grid.  The 
vertical  axis  gives  the  dynamic  residual.  The  MG  procedure 
with  three  levels  ( M  -  3)  shows  a  much  better  convergence 
rate  than  the  single  grid  solution  for  the  same  problem.  The 
convergence  factor  per  relaxation  iteration  (see  Ref.  10  for 
definition)  for  M  -  3  is  ji  -  0.60,  whereas  for  the  single  grid 
case  ( Af  -  1)  (i-  0.97.  The  linearized  PNS  equations  were 
also  solved  without  the  streamwise  correction  of  Refs.  3  and  4.  • 
The  MG  convergence  factor  is  much  worse  (/i  -  0.79)  Upon 
increasing  the  number  of  grid  points,  the  unmodified  equa¬ 
tion’s  convergence  deteriorates.  As  one  can  expect,  the  cor¬ 
rected  equauons  and  the  solution  of  the  equivalent  single 
Poisson  equation  for  the  pressure  exhibit  very  similar  conver¬ 
gence  properties. 

Nonlinear  Cases 

A  series  of  flow  problems  were  solved  with  the  nonlinear 
PNS  equations.  Several  test  cases  were  run  in  a  Cartesian 
coordinate  system  with  possible  clustering  of  mesh  points  by 
one-dimensional  stretching  functions  in  each  direction.  Among 
them  we  shall  mention  the  flow  over  an  infinite,  a  semi-infinite, 
and  a  finite  flat  plate,  trailing  edge  flow  and  an  entrance  flow 
between  two  flat  plates.  Two  cases  were  run  with  curvilinear 
orthogonal  coordinate  systems:  the  flow  along  an  axisymmet- 
ric  cylinder  and  the  flow  over  a  prolate  spheroid  at  zero  angle 
of  attack  In  the  following  sections  some  of  the  cases  will  be 
detailed 

Semi  Infinite  Flat  Plate 

In  this  case  the  flow  is  computed  starting  from  the  leading 
edge,  where  a  uniform  velocity  U  —  1,  F  -  0  is  given.  The 
downstream  boundary  was  set  at  Z  -  { **  4  and  zero  pressure 
gradient  was  specified  there  as  dP/dZ  -  0.  On  the  outer 
boundary,  U  -  l  and  P-0  The  no-slip  and  no-injection 
conditions  were  used  at  the  plate. 

The  second-order  accuracy  convergence  of  the  finite- 
difference  equations  is  demonstrated  in  Fig.  5.  The  pressure  at 
a  fixed  point  and  the  normalized  skin  friction  C)  —  CF<J Re  ■  Z 
(which  is  proportional  to  the  main  velocity  gradient)  at 
/  -  4  are  plotted  against  H2;  H  is  proportional  to  the  mesh 

interval. 

The  convergence  history  on  a  finest  grid  consisting  of 
65  X  65  points  is  shown  in  Fig.  6  for  Re  -  10J.  The  outer 
boundary  is  placed  at  Tro„  -  1  where  Y  is  the  normal  coordi¬ 
nate  (  T  -  r|)  Increasing  the  number  of  levels  M  increases  the 
convergence  rate  until  M  —  4.  The  convergence  factor  per 
relaxation  (p  -  0.562)  is  close  to  the  theoretical  value  obtained 
for  the  solution  of  a  single  Poisson  equation  with  the  SLOR 
method  tp  -  0  547). 

For  Re  -  10*  and  Z  -  3,  the  pressure  distribution  as  a 
function  of  the  Blasius  similarity  coordinate  [  -  yyjRe/Z  is 
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Fig.  7  The  pressure  profile  al  Z  ■  3  (semi-infinite  flat  plate). 
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Fig.  S  Normalized  skin  friction  coefficient  and  wake  centerline  veloc¬ 
ity  development  (trailing  edge  flow). 


shown  in  Fig  7.  The  distance  from  the  leading  edge  is  large 
enough  to  satisfy  the  Blasius  assumptions;  the  skin  friction 
coefficient  deviates  no  more  than  3%  from  the  theoretical 
value.  The  computed  pressure  is  compared  with  the  “Blasius 
pressure”  that  can  be  obtained  from  the  normal  momentum 
equation 

4ReZ(  P  -  Pw)  —/)  +  2f/"  (12) 

Pw  is  the  pressure  on  the  plate.  /  is  the  solution  of  the  Blasius 
boundary-layer  equation  for  a  flat  plate  with  zero  pressure 
gradient,  f  and  f"  are  the  first  and  second  derivatives  with 
respect  to  {.  The  agreement  between  the  present  computations 
and  Eq.  (12)  is  good. 

Trailing  Edge  Flow 

The  Reynolds  number  was  fixed  at  Re  —  10s  and  a  grid  of 
65  x  73  points  was  used.  The  flat  plate  with  zero  thickness 
occupies  the  interval  [0,1}  along  the  Z(-  f)  axis.  We  assumed 
that  the  interaction  is  limited  to  the  interval  0.5  <  Z<  1.5. 
The  Blasius  solution  was  specified  at  (he  upstream  boundary 
(Z  - 0.5)  and  symmetry  boundary  conditions  were  given  be¬ 
hind  the  plate.  The  Blasius  solution  was  approximated  by  the 
fourth-order  Karman-Pohlhausen  velocity  profile. 

Figure  8  shows  the  normalized  skin  friction  coefficient 
(C,-  Cf/Re)  on  the  plate  (Z<  1)  and  the  wake  centerline 
velocity  (UC,Z>  1).  The  agreement  with  the  interacting 
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boundarv  layers  (IBL)  solution  and  with  the  triple  deck  (TD)13 
revolt  is  satisfactory. 


Fntru*ue  Flo*  Between  Two  Flat  Plates 

Uniform  inlet  conditions  were  specified  at  the  entrance 
It  —  1.1—0)  while  the  usual  no-slip  and  no-injection  condi¬ 
tions  were  given  at  the  plates.  At  the  downstream  boundary 
the  pressure  gradient  was  calculated  by  assuming  a  fully 
developed  flow 


tit  is  the  rate  of  mass  flow. 

Figure  9  shows  the  development  of  the  centerline  velocity 
for  two  Reynolds  numbers,  fie  -  20  and  200,  and  for  a  grid  of 
41  x  101  points.  The  results  are  compared  with  the  full 
Navier-Stokes  computations  of  Ref.  14.  Again  the  agreement 
is  satisfactory,  even  at  the  entrance  where  the  omission  of  the 
streamwise  diffusion  terms  may  be  questionable. 

Prolate  Spheroid  at  Zero  Angle  of  Attack 

This  test  case  is  more  stringent  because  a  curvilinear  or¬ 
thogonal  prolate  spheroid  coordinate  system  is  used  and  the 
flowfield  is  more  complex  than  in  the  previous  cases.15  A 
nonzero  pressure  gradient  exists  at  the  outer  boundary: 
favorable  at  the  front  half  and  adverse  in  the  rear  half.  The 
flowfield  eventually  separates,  with  reversed  flow  near  the 
spheroid  Several  numerical  solutions  of  the  present  problem 
exist,15  ”  but  all  of  them  use  the  boundary-layer  approx¬ 
imation. 

The  flowfield  was  computed  for  a  region  between  the  nose 
and  the  rear  part  of  the  spheroid.  The  analytically  known 
potential  flow  is  specified  at  the  outer  boundary  (U  and  P) 
and  at  the  downstream  boundary  ( dP/di).  At  the  spheroid 
the  no-slip  and  no-injection  conditions  is  given.  At  the  up- 
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Fig-  9  The  development  ot  the  centerline  velocity  (entrance  flow). 
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Fig.  10  The  distribution  ot  the  skin  friction  coefllcient  tk>n|  the 
spheroid  of  thickness  4:1. 


stream  boundary  the  two  velocity  components  are  specified. 
Usually  they  can  be  computed  from  a  boundary-layer  code.  In 
the  present  work,  the  approximate  main  velocity  was  com¬ 
puted  from  a  Karman-Pohlhausen  profile  with  a  specified  skin 
friction  coefficient  or  displacement  thickness.  The  specified 
quantity  was  taken  from  existing  boundary-layer  solutions  of 
the  same  case.  The  normal  component  of  the  velocity  V  was 
determined  from  the  potential  solution.  Numerical  experi¬ 
ments  show  that  for  reasonable  upstream  conditions  the  solu¬ 
tion  is  independent  of  the  upstream  conditions,  apart  from  the 
first  few  marching  steps. 

The  downstream  boundary  was  set  before  the  separated 
zone,  although  the  PNS  equations  are  valid  there.  The  compu¬ 
tation  of  fiowfields  with  reversed  main  flow  needs  modifica¬ 
tions  in  the  approximation  of  the  convection  terms. 

Two  thickness  ratios  of  the  axes  were  considered— 4 : 1  and 
6:1.  In  each  case  the  Reynolds  number  was  Re  -  106  and  a 
single  grid  consisting  of  17  X  65  points  was  used  (in  the  i)  and 
£  directions  respectively). 

Thickness  Ratio  4 :  / 

The  dependence  of  the  normalized  skin  friction  coefficient 
(Cf  -  Cf>fRe)  on  the  axial  distance  (Z  -  cos£)  is  shown  in 
Fig.  10.  It  compares  well  with  the  boundary-layer  solution  of 
Ref.  16  but  disagrees  with  Wang’s  solution.15 

The  point  of  separation  is  determined  where  Cf  vanishes.  In 
the  present  case,  the  separation  point  was  found  to  be  at 


Fig.  It  The  pressure  distribution  along  the  spheroid  of  thickness  6:1. 
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Fig.  12  The  distribution  of  the  skin  friction  coefllcient  along  the 
spheroid  ol  thickness  6: 1. 


'.V.WV 
<_  ■ 


AV-'A"'.’ 


A  .-v  , 


j’-V-.v’A’. 


.  I 


MAY  1987 


FLOWS  BY  A  MARCHING  MULTIGRID  NONLINEAR  METHOD 


647 


r  *  c 

0  o< 


-10  -0  5  0  0  0  5  1  0 

z 

Fig.  1.7  The  development  ol  the  integral  thicknesses  along  the  spheroid 
of  thickness  6:1. 


7.  -  0  67.  which  is  in  excellent  agreement  with  Refs.  16  and 
17.  Wang15  has  computed  the  separation  point  at  Z  -  0.84. 

Thickness  Ratio  6.  / 

The  pressure  distribution  on  the  spheroid,  as  is  found  in  the 
PNS  solution,  is  compared  with  the  potential  solution  in  Fig. 
11.  The  agreement  is  very  good  except  near  the  downstream 
boundary  It  demonstrates  the  validity  of  the  boundary-layer 
approximation  for  the  computation  of  unseparated  flow  re¬ 
gions  over  slender  bodies.  The  slight  disagreement  near  the 
downstream  may  be  attributed  to  the  improper  value  being 
given  for  the  downstream  pressure  gradient.  The  potential 
pressure  gradient,  as  was  chosen  in  this  case,  is  not  a  good 
choice  near  regions  of  separation. 

The  skin  fnction  coefficient  distribution  is  compared  with  a 
boundary-layer  calculation  of  Schoenauer17  in  Fig.  12.  The 
agreement  is  good  except  near  the  upstream  and  downstream 
boundaries.  The  disagreement  in  the  upstream  region  prob¬ 
ably  results  from  unsatisfactory  initial  conditions  for  the 
velocities  The  separation  point  is  found  to  be  closer  to  the 
nose  ( Z  *>  0  77)  than  in  Ref.  17  (Z-0.87),  presumably  be¬ 
cause  of  the  greater  adverse  pressure  gradient  computed  from 
the  PNS  solution.  For  this  more  slender  body,  the  separation 
point  is  further  downstream  than  in  the  previous  case  (with 
the  thickness  ratio  4 : 1). 

Figure  13  gives  the  displacement  (8*)  and  the  momentum 
(<*>,)  thicknesses  dependence  on  Z.  The  values  computed  in 
the  present  work  are  higher  than  in  Ref.  17.  The  difference  is 
more  pronounced  near  the  separation  point. 

VI.  Conclusion 

In  the  present  study,  an  iterative  marching  method  for  the 
solution  of  the  PNS  equations  is  extended  to  curvilinear 
orthogonal  coordinate  systems.  Stable  second-order-accurate 
timte-difference  equations  were  developed  using  central  ap¬ 
proximation  to  all  derivatives.  The  solution  algorithm  is  basi- 
eallv  equivalent  to  the  solution  of  a  single  elliptic  equation 
with  the  SLOR  method  The  efficient  two-dimensional  MG 
procedure  previously  developed  for  a  linearized  PNS  equation 
was  applied  to  the  solution  of  several  nonlinear  viscous  flows 


in  Cartesian  grids.  The  convergence  rates  were  not  adversely 
affected  by  the  nonlinearity. 
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